

% HHVR_EM1a_2_types.m is called by HHVR_2016_Master_code.m
% HHVR_EM1a_2_types.m calls  HHVR_EM1b_2_types.m 

  function  [loglik  ] = HHVR_EM1b_2_types(bbeta)  
  global DATA  X choices y_short zz individuals  
 
 
s=size(DATA,1)/choices;      
nume1=zeros(s,1);nume2=zeros(s,1);nume3=zeros(s,1) ;nume4=zeros(s,1) ;nume5=zeros(s,1) ; 
numu1=zeros(s,1);numu2=zeros(s,1);numu3=zeros(s,1) ;numu4=zeros(s,1) ;numu5=zeros(s,1) ; 
prob_e1=zeros(s,1);prob_e2=zeros(s,1);prob_e3=zeros(s,1); prob_e4=zeros(s,1);prob_e5=zeros(s,1); 
prob_u1=zeros(s,1);prob_u2=zeros(s,1);prob_u3=zeros(s,1); prob_u4=zeros(s,1);prob_u5=zeros(s,1);  
denu= zeros(s,1); dene= zeros(s,1);   
denu_x1= zeros(s,1); dene_x1= zeros(s,1);  

ze=zz(:,1); zu=zz(:,2); 
for i=1:(size(DATA,1))/choices                    % ie from 1 to number of individuals
       
     nume1 (i) =  exp(X(i*choices-4,1)*bbeta(5)    ) ;
     nume2 (i) =  exp(X(i*choices-3,1)*bbeta(5)  + 1*bbeta(1)  ) ;
     nume3 (i) =  exp(X(i*choices-2,1)*bbeta(5)  + 1*bbeta(2)  ) ;
     nume4 (i) =  exp(X(i*choices-1,1)*bbeta(5)  + 1*bbeta(3)  ) ;
     nume5 (i) =  exp(X(i*choices  ,1)*bbeta(5)  + 1*bbeta(4)  ) ;
     dene(i)   =  nume1 (i)  +  nume2 (i) + nume3 (i)+ nume4 (i)+ nume5 (i) ;
     dene_x1(i)   =  nume1(i)*X(i*choices-4,1)  +  nume2(i)*X(i*choices-3,1) + nume3(i)*X(i*choices-2,1)+ ... 
                     nume4(i)*X(i*choices-1,1) +   nume5(i)*X(i*choices  ,1 ) ;
     prob_e1(i)     =  nume1 (i) /  dene(i)  ;     prob_e2(i)     =  nume2 (i) /  dene(i)  ;
     prob_e3(i)     =  nume3 (i) /  dene(i)  ;     prob_e4(i)     =  nume4 (i) /  dene(i)  ;
     prob_e5(i)     =  nume5 (i) /  dene(i)  ;  
    
     numu1 (i) =  exp(X(i*choices-4,2)*bbeta(6)    ) ;
     numu2 (i) =  exp(X(i*choices-3,2)*bbeta(6)  + 1*bbeta(1)  ) ;
     numu3 (i) =  exp(X(i*choices-2,2)*bbeta(6)  + 1*bbeta(2)  ) ;
     numu4 (i) =  exp(X(i*choices-1,2)*bbeta(6)  + 1*bbeta(3)  ) ;
     numu5 (i) =  exp(X(i*choices  ,2)*bbeta(6)  + 1*bbeta(4)  ) ;
     denu(i)   =  numu1 (i)  +  numu2 (i) + numu3 (i)+ numu4 (i)+ numu5 (i) ;
     denu_x1(i)   =  numu1(i)*X(i*choices-4,2)  +  numu2(i)*X(i*choices-3,2) + numu3(i)*X(i*choices-2,2)+ ... 
                     numu4(i)*X(i*choices-1,2) +   numu5(i)*X(i*choices  ,2 ) ;
     prob_u1(i)     =  numu1 (i) /  denu(i)  ;     prob_u2(i)     =  numu2 (i) /  denu(i)  ;
     prob_u3(i)     =  numu3 (i) /  denu(i)  ;     prob_u4(i)     =  numu4 (i) /  denu(i)  ;
     prob_u5(i)     =  numu5 (i) /  denu(i)  ;  
 end
 
for i=1:(size(DATA ,1))/choices     
        
          if y_short(i)==1
          
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,3)- zu(i)*( prob_u2(i) )*X(i*choices,3)   ;
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,4)- zu(i)*( prob_u3(i) )*X(i*choices,4) ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,5)- zu(i)*( prob_u4(i) )*X(i*choices,5)   ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,6)- zu(i)*( prob_u5(i) )*X(i*choices,6)  ;
            
           score(i,5) =     ze(i,1)*(X(i*choices-4,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-4,2)- (denu_x1(i) / denu(i)));  
            
          elseif y_short(i)==2
  
           score(i,1) =     ze(i)*(1-prob_e2(i))*X(i*choices,3) + zu(i)*(1-prob_u2(i))*X(i*choices,3)  ;
           
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,4)- zu(i)*( prob_u3(i) )*X(i*choices,4) ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,5)- zu(i)*( prob_u4(i) )*X(i*choices,5)  ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,6)- zu(i)*( prob_u5(i) )*X(i*choices,6)  ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-3,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-3,2)- (denu_x1(i) / denu(i)));   
           
           elseif y_short(i)==3
  
           score(i,2) =     ze(i)*(1-prob_e3(i))*X(i*choices,4) + zu(i)*(1-prob_u3(i))*X(i*choices,4)  ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,3)- zu(i)*( prob_u2(i) )*X(i*choices,3)   ; 
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,5)- zu(i)*( prob_u4(i) )*X(i*choices,5)   ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,6)- zu(i)*( prob_u5(i) )*X(i*choices,6)  ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-2,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-2,2)- (denu_x1(i) / denu(i)));   
           elseif y_short(i)==4
  
           score(i,3) =     ze(i)*(1-prob_e4(i))*X(i*choices,5) + zu(i)*(1-prob_u4(i))*X(i*choices,5)   ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,3)- zu(i)*( prob_u2(i) )*X(i*choices,3)  ; 
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,4)- zu(i)*( prob_u3(i) )*X(i*choices,4)  ;
           score(i,4) =   - ze(i)*( prob_e5(i) )*X(i*choices,6)- zu(i)*( prob_u5(i) )*X(i*choices,6)  ;
             
           score(i,5) =     ze(i,1)*(X(i*choices-1,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices-1,2)- (denu_x1(i) / denu(i))); 

           elseif y_short(i)==5
  
           score(i,4) =     ze(i)*(1-prob_e5(i))*X(i*choices,6) + zu(i)*(1-prob_u5(i))*X(i*choices,6)  ;
           
           score(i,1) =   - ze(i)*( prob_e2(i) )*X(i*choices,3)- zu(i)*( prob_u2(i) )*X(i*choices,3)   ; 
           score(i,2) =   - ze(i)*( prob_e3(i) )*X(i*choices,4)- zu(i)*( prob_u3(i) )*X(i*choices,4)  ;
           score(i,3) =   - ze(i)*( prob_e4(i) )*X(i*choices,5)- zu(i)*( prob_u4(i) )*X(i*choices,5)   ;
             
           score(i,5) =     ze(i,1)*(X(i*choices ,1)- (dene_x1(i) / dene(i)));
           score(i,6) =     zu(i,1)*(X(i*choices ,2)- (denu_x1(i) / denu(i)));  
          else
              disp 'not good'
              pause
          end
     
end
  
sscore = sum(score) ;
loglik = [sscore ];
  
  
  
  